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Hubbard ladders are an important stepping stone to the physics of the two-dimensional Hub¬ 
bard model. While many of their properties are accessible to numerical and analytical techniques, 
the question of whether weakly hole-doped Hubbard ladders are dominated by superconducting or 
charge-density-wave correlations has so far eluded a dehnitive answer. In particular, previous numer¬ 
ical simulations of Hubbard ladders have seen a much faster decay of superconducting correlations 
than expected based on analytical arguments. We revisit this question using a state-of-the-art imple¬ 
mentation of the density matrix renormalization group algorithm that allows us to simulate larger 
system sizes with higher accuracy than before. Performing careful extrapolations of the results, 
we obtain improved estimates for the Luttinger liquid parameter and the correlation functions at 
long distances. Our results confirm that, as suggested by analytical considerations, superconducting 
correlations become dominant in the limit of very small doping. 

PACS numbers: 71.27.+a, 74.20.Rp, 74.72.Gh, 02.70.-c 


I. INTRODUCTION 

The question of whether electrons in two dimensions 
can exhibit superconductivity mediated by repulsive in¬ 
teractions, which is motivated by the discovery of high- 
temperature superconductors, has become one of the cen¬ 
tral questions of condensed matter theory. However, the 
numerical study of even the simplest models, such as the 
Hubbard or t-J model, are made difficult by a multitude 
of competing low-energy phases that these models ex¬ 
hibit. This is particularly the case in the regime of weak 
doping away from half filling, which is most relevant for 
the phase diagram of cuprate superconductors. Numer¬ 
ical efforts reviewed in Ref. 1 as well as more recent re¬ 
sults in Refs. 2 and 3 have shown a close competition of 
striped antiferromagnetic phases, d-wave superconduct¬ 
ing phases, and other more exotic phases such as a pseu¬ 
dogap phase where hole quasi-particles play the role of 
the mobile carriers. 

Faced with these challenges, quasi-one-dimensional 
systems such as ladders have appeared as an easier start¬ 
ing point to investigate the properties of these models, as 
they are amenable to a broader range of numerical and 
analytical methods. These approaches view the system 
as essentially one-dimensional with additional degrees of 
freedom that allow the two-dimensional characteristics 
to emerge. Crucially, this has allowed treatment using 
the density matrix renormalization group (DMRG),^’® 
which allows accurate simulations of extended quasi- 
one-dimensional systems and has successfully illuminated 
many properties of ladder systems.® 

Numerical work on t-J and Hubbard ladders^”® as well 
as analytical work on the weak-interaction limit by Ba- 
lents and Fisher^® has shown that in a wide parameter 
regime, weakly doped ladders fall into the Luther-Emery 
universality class,which has a gapped spin mode and 
a single gapless charge mode. This phase is a possible pre¬ 


cursor phase to two different ordered phases in the two- 
dimensional limit, a superconducting phase (SC) and a 
charge-density wave (CDW) phase. To distinguish these 
phases, it is crucial to compute whether the ladder system 
is dominated by density-density or superconducting cor¬ 
relations. Within the Luther-Emery universality class, 
these both decay with a power-law whose exponents are 
determined by a single dimensionless parameter Kp. 

Previous DMRG calculations of the correlation 
functions^®’^^ have observed power-law decay of the cor¬ 
relation functions, but have found a surprisingly fast de¬ 
cay of the superconducting correlations inconsistent with 
dominant superconducting correlations. This decay was 
found to be inconsistent with calculations in the weak- 
doping limit,and also violates certain identities of 
the Luther-Emery universality class. Here, we revisit the 
calculation of these exponents with a focus on extracting 
the correct behavior of the pair correlation function. We 
exploit the advances in DMRG methods, in particular on 
the correct extrapolation of physical quantities, and in¬ 
creases in computational power since the work of Refs. 13 
and 14. Using a high-performance DMRG code,^^ we are 
able to target longer systems with much improved ac¬ 
curacy to obtain reliable correlation exponents for the 
two-leg Hubbard ladder, which settle the disagreement 
between the numerical calculations and the theoretical 
expectations. To achieve this, we carefully analyze the 
effects of the hnite system size and the DMRG truncation 
on correlation functions. 

The paper is organized as follows. In Section H we 
present the Hubbard model. In Sections HI we briefly 
introduce the DMRG method and how correlation ob¬ 
servables are extrapolated to the thermodynamic limit. 
Sections IV and V are devoted to the discussion of our 
results: first the finite size analysis of density oscilla¬ 
tions, then the comparison of the Luttinger liquid expo¬ 
nent with the pair and density correlation functions. In 
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r = \i- j\ 


Figure 1. (color online) Spatial decay of the pair correlation 
function (blue markers) for a 2 x 32 Hubbard model with 
U/t = 8 and average filling n = 0.875 similar to the results 
obtained in Ref. 14. The solid, dotted and dashed lines are 
reference power-law decays with exponents /r = —1/2,—1,-2, 
respectively. In this paper we estimate an exponent /r « — 1 
(see exponents in Table I). 


Section VI we present our conclusions. 


II. THE HUBBARD LADDER 


We consider the Hubbard model on a two-leg ladder 
described by the Hamiltonian 


H — t ^ + H-C. 

•i,A,cr 

i,(T 

V U ^ ^ 

2, A 


( 1 ) 


where the index i runs along two coupled chains of length 
L and A = 1, 2 identifies the two chains respectively. The 
operator ^ creates a fermion at site i on chain A 

with spin tr S {t,i} and h(i^x),a = 

The phases of this model can conveniently be labeled 
by the number of gapless spin and charge modes, with 
up to two gapless modes possible in each sector. Much 
attention has been focused on the phase with one gap¬ 
less charge mode, but a gap in the spin sector (labeled 
CISO in Ref. 10) for its relevance for both superconduct¬ 
ing (SC) and charge-density wave (CDW) phases in the 
two-dimensional limit. This phase is found in a wide pa¬ 
rameter range of repulsive U, t± < 2t, and hole-doping 
with filling n < 1 (in units where one fermion per site 
corresponds to n = 1). In this paper we will focus on the 
isotropic hopping case t± = t with interaction U/t = 8, 


for which the spin gap has previously been reportedto 
show a maximum. We investigate different values of the 
average filling n while keeping the total magnetization 
fixed at zero. 

We define the local rung density operator as = 
a ^(bA),o- ^n.d its expectation value as 

~ ' (b-(2,A),g)- (2) 

A,fT 

Its density correlation function is 

N{i,j) = {n.rij) - (3) 

and the d-wave pair correlation function takes the form 

iA(*,j) = (AlA,), (4) 

where AJ = “ 4.i)d creates a singlet 

on rung i. 

In the Luther-Emery phase, the spatial decay of the 
density-density correlation function N{r) and the pair 
correlation function D{r) at large distance r are domi¬ 
nated by a power-law parametrized by the non-universal 
parameter Kp-. 

N{r) oc r~'' with v = Kp, (5) 

D{r) oc r~^ with fi = l/Kp. (6) 

Because of the relation • /i = 1, one has that for Kp > 1 
the system is dominated by the d-wave pair correlations, 
whereas for Kp < 1 one observes dominant charge den¬ 
sity wave correlations. The Luttinger liquid parameter 
Kp must in general be determined numerically. In the 
limit n —> 1 and in the strong-coupling limit of the t-J 
model, one can construct an effective bosonic model for 
hole pairs in the Hubbard model, -^^hich yields a uni¬ 
versal power-law decay of the pair correlation function 
Dir) oc ^l\fr, and thus Kp = 2. Previous DMRG cal¬ 
culations, whose results we reproduce in Fig. 1, found 
the decay of the pair-correlation function to be much 
faster than Dir) cx 1/v/r for weak doping, and a compar¬ 
ison of the decay of pair- and density-density correlations 
showed a violation of the identity v ■ fi = 1. 


HI. SIMULATION METHOD 

A. The density matrix renormalization group 
algorithm 

We tackle the model using an implementation of 
the density matrix renormalization group (DMRG) 
method^’”’ in a formalism of matrix-product states 
(MPS)^®d9 available as part of the ALPS project. 

For recent reviews of these methods, see Refs. 6 and 21. 

The DMRG method can be understood as a variational 
optimization over Matrix Product State (MPS) wave- 
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Figure 2. (color online) Spatial decay of the pair correlation 
function D{r) on ladders with L = 128 and n = 0.875 as a 
function of distance r = \i — j\ for several choices of i: fixing 
i = 18 (orange line) and i = 20 (green line), and averaging 
11 pairs {i,j) at distance r around the middle according to 
Eq. (7) (blue line). 


functions, which are a class of one-dimensional ansatz 
states that can be systematically refined by increasing 
the so-called bond dimension M. In a standard approach, 
the variational optimization proceeds by iteratively im¬ 
proving the wavefunction on pairs of sites. In each op¬ 
timization step, a truncation occurs, and the sum of the 
discarded components of the wave function, called trun¬ 
cated weight e, is stored for later evaluations (see Sec¬ 
tion III C). 


1.95 


1.90 

1.85 

e 

1.80 

1.75 

1.70 


0 20 40 60 80 100 120 

rung i 

0—0 M = 800 V—V M = 2000 <—< M = 3600 

o—□ M = 1200 ^ M = 2800 o—o M = 4000 

M = 1600 0—0 M = 3200 □—a M = oo 



Figure 3. (color online) Local density profile for ladders with 
L = 128 and n = 0.875. for several bond dimensions M, 
the inset showing details in the center of the system. The 
uncertainty due to systematic errors in the extrapolated curve 
M = oo is shown as a shaded region around the estimate; its 
size is often smaller than the symbols. 


the middle of the system 


C{r) 
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( 7 ) 


While the MPS ansatz is exact for M —>■ oo, a finite 
value of M restricts the amount of entanglement that can 
be captured by the wave function. It has been shown^^^^^ 
that this is an efficient representation of the ground states 
of one-dimensional, gapped, local Hamiltonians. For gap¬ 
less systems with a dynamical critical exponent of z = 1, 
one finds that only a polynomially growing bond dimen¬ 
sion^®’^® is generally required to accurately describe local 
properties. When coupling chains to ladders, M has to 
increase exponentially with the width of the ladders, as 
the entanglement entropy grows linearly with the width. 

DMRG is most efficiently performed with open bound¬ 
ary conditions, which causes complications since this 
breaks translational invariance. Local quantities (for ex¬ 
ample density or magnetization) will differ from site to 
site and correlation functions depend not only on the dis¬ 
tance r = \i — j\ between two sites but on both sites i 
and j. The latter effect can be observed in Fig. 2, where 
we plot D{r) at various distances between one fixed site 
i and an other site j = i + r. To reduce boundary effects 
for a correlation function C{r), we evaluate the value at 
distance r by averaging over 11 pairs r = \i — j\ around 


The solid dark line in Fig. 2 shows the impact of averag¬ 
ing for reducing the oscillations induced by the bound¬ 
aries. 

In the simulations shown here, we use a two-site up¬ 
date algorithm and generally perform between 20 and 30 
sweeps until energy and local observables converge for 
a given bond dimension M. For all model parameters, 
we perform independent simulations with several bond 
dimensions up to a maximum oi M = 4800. 


B. Finite size and finite entanglement scaling 

To obtain reliable long-range correlation functions, 
DMRG results need to be extrapolated both in system 
size L and bond dimension M. In Figs. 3 and 4 we show 
the local density and the pair correlation functions for 
various values of M and L. It is apparent in Fig. 4 that 
calculations with insufficient system size or bond dimen¬ 
sion may lead to underestimating the strength of the cor¬ 
relations. 

To understand the interplay of bond dimension and 
system size, it is crucial to note that a matrix-product 
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Figure 4. (color online) Spatial decay of the pair correlation 
function D{r) for ladders with n = 0.875 and several bond 
dimensions M and system sizes L. Top panel: Results for 
L = 128. Bottom panel: Results are extrapolated to M = 
oo. The shaded region around the curve show the confidence 
range of the systematic error. 


state always exhibits exponentially decaying correlations 
at large enough distances,'® but may reproduce power- 
law decay at short distances. The scale on which the 
correlations cross over from power-law to exponential be¬ 
havior, is dictated by the bond dimension M. At 
the same time, the finite size of the system will introduce 
some length scale on which boundary and finite-size 
effects become significant, and correlations are no longer 
representative of the thermodynamic limit. In interpret¬ 
ing results obtained with DMRG for finite M and L, it 
is important to distinguish different regimes depending 
on whether deviations from the asymptotic behavior are 
dominated by or by Following Ref. 27, the regime 
of small bond dimension M, where ^ is referred 
to as finite entanglement scaling (FES) regime, in which 
the system does not feel the presence of the boundaries 
because the correlation length induced by the finite bond 
dimension is short ranged compared to the system size. 
In the other limit, where exhausts and correla¬ 
tions can in principle span the whole system, the finite 
size scaling (FSS) regime is reached. In an intermedi¬ 
ate regime, where and are comparable, a two- 
parameter scaling may be necessary. 

To illustrate these regimes, the top panel of Fig. 4 
shows D{r) for several bond dimensions M at L = 128. 


We can clearly distinguish how correlations are cut off at 
a certain length scale depending on the bond dimen¬ 
sion; as the bond dimension M is increased, we can con¬ 
sider the correlations converged for the given system size 
over an increasing range of distances. The bottom panel 
shows results that have been extrapolated in M for dif¬ 
ferent system sizes L, and thus suffer only from finite-size 
corrections. By comparing the two panels, we see that for 
the range of system sizes considered here, ~ 50 and 
thus the results in the upper panel suffer primarily from 
finite-entanglement corrections for r < 50. For the bond 
dimensions we can attain in practice, corrections due to 
set in at shorter distances, as seen in the upper panel. 

In light of these considerations, we avoid having to 
perform a two-parameter scaling and focus mostly on the 
more relevant corrections due to finite entanglement. We 
thus perform careful extrapolations in the bond dimen¬ 
sion M for a fixed, given system size, achieving the finite- 
size scaling limit ^ and then compare results ob¬ 
tained for different system sizes to assess the reliability. 
Most of the results in the remainder of this paper are 
obtained for L — 128. 


C. Extrapolation to infinite bond dimension 

When extrapolating observables to infinite bond di¬ 
mension M = oo one choice is to extrapolate in 1/M. 
However, a more controlled extrapolation may be possi¬ 
ble by extrapolating in the variance of the energy 

Var[i7] = (H2) _ {H)\ (8) 

which vanishes for an eigenstate of H. Similarly, the 
truncated weight e will vanish when the bond dimen¬ 
sion M is large enough to faithfully represent the wave 
function. For a more reliable data analysis we com¬ 
pare extrapolations in both 1/M, truncated weight e and 
Var[IT]. Unless noted otherwise, we will plot the average 
of the three extrapolations together with a confidence in¬ 
terval given by the minimum and maximum extrapolated 
values as an estimate of the systematic error. When fur¬ 
ther analysis is performed, e.g. for the determination of 
ATp, the analysis is performed for both extrapolations to 
obtain error estimates on the results. 


1. Extrapolating the ground state energy 

For the ground state energy, the deviation of {H) from 
the ground state energy Eq is known^® to depend linearly 
on Var[M], which provides for very accurate extrapola¬ 
tions. To demonstrate this dependence we write the state 
IV') = IV'o) + |<5) obtained by DMRG as the sum of the 
true ground state |V’o) with energy Eq and an error term 
|5) with (5|(5) = (5^. Both |V’) and |V’o) are supposed to 
be normalized. The energy of this state then is (with 
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Figure 5. (color online) Extrapolation of the ground state 
energy density as a function of the inverse bond dimension 
1/M (left panel), the truncated weight e (right panel, full 
circles) and as a function of the energy variance Var[i^] (right 
panel, empty circles). Results are for L = 128 and n = 0.875. 


5 = {tpo\S) + ((5| V^o)) 

{H) = {i;\H\^)=Eo{l + S) + 0{S^). (9) 

Similarly, the expectation value of is 

{H^)=El{\ + 5)^0{6‘^). ( 10 ) 

Combining these results we obtain for the energy vari¬ 
ance (8): 


{H^)-{Hf = El5 + 0{5^), ( 11 ) 

which can be used to derive a linear dependence of the 
expectation value of the ground state energy on the vari¬ 
ance: 


{H) =Eo + aVar[77] -h 0{5^), (12) 

where a is a non-universal pre-factor independent of the 
DMRG error 5. This linear dependence can be seen in 
Fig. 5, and provides a more reliable extrapolation than 
by extrapolating naively in 1/M. 




1000/M Var[.tf] or e x 10^ 


2. Extrapolating other observables 

This strategy cannot be generalized to generic observ¬ 
ables that do not commute with the Hamiltonian. We 
thus perform linear regressions with quadratic polyno¬ 
mials in both the energy variance, the truncated weight 
and 1/M, using results for the six largest values of M. 
Spatially dependent quantities, such as the local density 
and the correlation functions, are independently extrap¬ 
olated at each point. As a consistency test sum rules are 
checked, e.g. the sum of all local densities should give 
the total particle number. 

Figure 6 shows the extrapolations for N{r) and 
D{r) at various positions using both extrapolation 


Figure 6. (color online) Extrapolation of observables as a 
function of the inverse bond dimension 1/M (left panels), the 
truncated weight e (right panels, full symbols) and as a func¬ 
tion of the energy variance Var)!?] (right panels, empty sym¬ 
bols). Shown are representative examples of (a) the local rung 
density Ui, (b) the density correlation function N{r), and (c) 
the pair correlation function D{r). Results are for L = 128 
and n = 0.875. 


schemes. One notices at first that the two approaches 
agree reasonably with the disagreement limited to a few 
percent. For local observables, such as the local density 
in panel (a), the difference can be very small, and de¬ 
viations are generally found to be smaller towards the 
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edge of the system. Correlation functions like N{r) and 
D{r) are harder to extrapolate as r increases and the 
relative differences in extrapolated values grow. Accu¬ 
rately describing the behavior of correlation functions at 
longer distances requires an increasing bond dimension 
M. Higher order terms in the extrapolation thus become 
more important at constant M and extrapolations are 
harder. 


IV. DETERMINING Kp FROM DENSITY 
OSCILLATIONS 


The most reliable estimation of the correlation expo¬ 
nent Kp in DMRG calculations is based on density os¬ 
cillations (Friedel oscillations) induced by the boundaries 
of the open systems commonly studied in DMRG.^® This 
method has been successfully applied to t-J ladders,^® 
quantum wires'^'^ and supersymmetric fermions. 

Friedel oscillations observed in the local density profile 
take the form 


n(x) = {n{x)) 


cos{2TTNhx/Les + (t>i) 

-R- /9 “f 

[Leff sin(7ra:/Leff + <1)2)] ” 


( 13 ) 

where A is a non-universal amplitude, (f>i and (f >2 phase 
shifts, no the background density, Nh is the number of 
holes in the system and Leff ~ L is an effective length. 
To derive this expression. Ref. 29 considered the slow¬ 
est decaying component of the density-density correlation 
function in the Luther-Emery model, finite size effects 
are then introduced via a standard conformal transfor¬ 
mation. 

The need for an effective length can be understood 
as an effect of the finite extent of the hole pairs in the 
ladder. We find that an effective length Leg = L — 2 best 
describes our results. For our choice of doping the ladders 
with multiples of four holes, there is a density maximum 
in the center of the ladder and we thus need to fix the 
phase shifts to (/>! = — 7r7V/iL/Leff and ^2 = f (1—-^^/Aeff)- 

The exponent Kp can be obtained from finite size scal¬ 
ing at constant density p. In particular, the density os¬ 
cillation in the center of the system follow 


6n{L) = n{L/2) - no ^ (14) 


Unfortunately, extracting the oscillation amplitude is 
not as straight-forward. The first problem is that the 
background density ng is not simply the mean density 
n, but instead depends on L. This can be seen by in¬ 
tegrating Eq. (13) and equating it to the total particle 
number. Numerically, this deviation can be observed in 
the top panel of Fig. 7. Secondly, for an even number of 
rungs, the finite lattice spacing limits the spatial resolu¬ 
tion, hence often it is not possible to obtain the density 
exactly at x = L/2, where the oscillating factors are triv¬ 
ial. This is not the case for an odd number of rungs as 
one can see in Fig. 7b. From the lower panel of Fig. 7 


(a) 



(b) 




Figure 7. (color online) (a) Fit of the density profile (solid 
line) for a system of length L = 128 with filling n = 0.875 
compared to the raw density (markers). The red dotted line 
is the density offset no obtained from the least-square fit, as 
a reference, the dashed green line shows the average filling n. 
The fit is restricted to the shaded region to avoid the divergent 
boundaries, (b) Same as (a) for a system with an additional 
rung but the same number of holes, i.e. with N-f = N_i = 113 
(two more particles and the same number of holes), (c) Finite 
size scaling of the oscillation amplitude 5n{L) as a function of 
the system size L in double logarithmic plot for several filling 
parameters. 
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Table I. Kp as extracted from the fit in Fig. 7. 
Var[77] e 1/M 


n 

Kp 


Kp 


Kp 


0.875 

0.99 

0.9999 

0.92 

0.9999 

1.17 

0.9956 

0.9375 

1.54 

0.9919 

1.54 

0.9916 

1.66 

0.9973 


we note that the latter problem has a minor impact on 
the final result, therefore we continue the analysis with 
ladders of even length. 

To avoid the first issue we perform a non-linear least 
square fit of the local densities in the middle of our system 
to Eq. (13). The parameters Teff, (/i and 4)2 are fixed as 
discussed above and A, Kp and no are used as fit param¬ 
eters. The obtained fit is then used to compute dn{L), 
and Kp is extracted in another fit to Sn{L) ~ 

This approach works very well, as one can see from the 
illustrative examples in Fig. 7. Here, we also show the 
difference between the size-dependent background den¬ 
sity no and the mean density n. 

Other approaches for extracting the oscillation ampli¬ 
tude provided no reliable results. One of the failed at¬ 
tempts was to fix no from a linear interpolation of the 
two closest points to the nodes in the oscillations. We 
also tried obtaining n(L/2), by accounting for the sin 
and cos terms in the finite size scaling, but this approach 
was unstable because of numerical errors in obtaining the 
wavelength of the oscillations. 

Our results for the exponent Kp obtained from the fit 
in Fig. 7c) are summarized in Tab. I. We see that Kp in¬ 
creases with filling and it is consistent with reaching the 
limit Kp = 2 for n = 1. The goodness of the linear regres¬ 
sion E? is reported to be always larger than 99%, which 
supports the expected decay of the oscillations. Results 
for more dilute systems are discussed in Appendix A. 

V. CORRELATION FUNCTIONS 

Our results at low doping are consistent with the ex¬ 
pected dominant superconducting correlations. The cor¬ 
relation exponent obeys Kp > 1, and therefore super¬ 
conducting pair correlations dominate and decay slower 
than 1/r. This is in contrast to Fig. 4 of Ref. 14, where 
the pair correlations decay roughly as r~^-^ at n = 0.875 
and n = 0.9375. 

This puzzling discrepancy can be resolved by noting 
that it is hard to faithfully obtain the correlation expo¬ 
nent from the spatial decay of correlation functions. In 
particular, as we argued in Sec. Ill, large system sizes L 
and large bond dimensions M are needed to obtain reli¬ 
able results for correlation functions. Both finite M and 
finite L suppress long range correlations. While the sys¬ 
tem sizes and bond dimensions of Ref. 14 are sufficient to 
converge local quantities, they are insufficient for corre¬ 
lation functions, and in particular the pair correlations. 


(a) 



r=\i-j 


(b) 



r=\i-j 


Figure 8. (color online) Decay of the density-density correla¬ 
tion function N{r) for many system system sizes. The dashed 
line is a power-law decay obtained with an exponent —Kp and 
Kp = 1.54 from the Var[i7] and e extrapolations in Table I; 
the vertical offset is chosen to show the agreement with the 
long-distance behavior of the correlation function. Panels a) 
and b) refer to two different average fillings n = 0.875 and 
n = 0.9375, respectively. 


Since the determination of Kp from local densities via 
Friedel oscillations is expected to be more reliable than 
from the decay of the correlation functions, we here use 
the values obtained in the previous section and show that 
the correlation functions - after proper extrapolation in 
M and L - are consistent with these values. Figures 8 
and 9 show N{r) and D{r) for various system sizes. We 
notice three regimes: At short distances we find a non- 
universal regime of fast-decaying correlation functions. 
At the longest distances finite-size effects become relevant 
and the correlation functions are strongly suppressed. In 
between we find a region where the spatial decay of the 
correlation functions is indeed consistent with a power 
law. For L = 192 and extrapolating M —>■ oo we find 
good agreement in the range 10 ;< r < 90 with the ex¬ 
pected behavior based on the values of Kp obtained by 
















(a) 


q 




faster than expected. The first cause is the need for 
very long system sizes, as finite sizes tend to strongly 
suppress pair correlation functions when the distance be¬ 
comes comparable to the system size. More importantly, 
a careful extrapolation in the bond dimension to M —>■ oo 
is necessary and has to be performed separately for each 
distance. Increasingly larger bond dimensions M are 
needed to obtain converged results for longer distances 
r, substantially larger than necessary to converge local 
quantities. 

We devote particular attention to the extrapolation 
techniques. All data used for our work and fitting and 
extrapolation workflows are available as Supplementary 
Material*^ to allow the reader to reproduce our results 
and modify the details of the extrapolation and fit ap¬ 
proaches and see how they affect the final results. 

While we here use a standard finite-size DMRG ap¬ 
proach with open boundary conditions, recently proposed 
techniques such as the sine-square deformation*^, grand- 
canonical or infinite-size DMRG'^® could also 

be applied to this problem. 

Achieving a good understanding of the two-leg lad¬ 
der and the effects of finite entanglement scaling and 
finite size scaling of DMRG observables, and obtaining 
reliable results for the pair correlation function in this 
simple model is an important milestone to improving 
the reliability of numerical simulations for larger, two- 
dimensional systems. 
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Figure 9. (color online) Decay of the pair correlation function 
D(r) for many system system sizes. The dashed line is a 
power-law decay obtained with an exponent — IjKp and Kp = 
1.54 from the Var[//] and e extrapolations in Table I; the 
vertical offset is chosen to show the agreement with the long¬ 
distance behavior of the correlation function. Panels a) and 
b) refer to two different average fillings n = 0.875 and n = 
0.9375, respectively. 

the extrapolations in Va,r[H] and e in Tab. I. 

Note that where the pairfield correlations always re¬ 
main positive, the density correlations oscillate across 
zero. This eventually leads to the spikes when N{r) is 
about to change sign in Fig. 8, where we plot its absolute 
value on a double-logarithmic axis. 

VI. CONCLUSIONS 

In this paper we settle the long-standing disagreement 
between the analytically predicted behavior of the pair 
correlation functions in weakly doped Hubbard ladders 
and results of DMRG calculations. We illustrate the two 
main causes of the discrepancy in previous results, which 
had indicated that the pair correlation function decays 
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Appendix A: Weak hole-doping results 

For very dilute systems such as n = 0.96875 we find 
that the finite size analysis described in Section IV be¬ 
comes less reliable, as DMRG convergence and extrap¬ 
olation become more challenging. Since we evaluate an 
even number of hole pairs, only three system sizes are 
available, L = 64, L = 128 and L = 192, correspond¬ 
ing to only 2, 4 and 6 hole pairs, respectively. Longer 
system sizes would be needed to perform a rigorous scal¬ 
ing analysis. Furthermore, distributing six hole pairs in 
such a long system is a very slow process and often leads 
to convergence problems as has been observed for the t-J 
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Table II. Kp as extracted from the fit in Fig. 10, where only 
the two largest system sizes L = 128 and L — 192 are consid¬ 
ered. 


n 

Var[i7] 

e 

1/M 

0.96875 

1.87 

1.85 

2.39 


model.To improve convergence speed one could employ 
multigrid techniques.'^® 

Here we present the results obtained at average filling 
n = 0.96875. Figure 10 shows that finite size scaling anal¬ 
ysis, whose exponents are reported in Table II. Compar¬ 
isons for the spatial decay of correlation functions with 
the exponent Kp expected from the Friedel oscillations 
are shown in Figure 11 and 12. 

Note that in the scaling analisys of the Friedel os¬ 
cillations we consider only the two largest system sizes 
L = 128 and L = 192, because the first system size 
L = 64 contains only two hole pairs, hence finite size ef¬ 
fects are expected to have a dominant contribution. The 
value of Kp obtained from the fit (see Table II) is again 
compatible with the expected limit for n = 1, and it is 


compatible with the decay of correlation functions (see 
Figure 11 and 12). 



Figure 10. (color online) Finite size scaling of the oscillation 
amptitude Sn{L) as a function of the system size L in double 
logarithmic plot for several an average filling n = 0.96875. 
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r = \i- j\ 


Figure 11. (color online) Decay of the density-density correla¬ 
tion function N (r) for many system system sizes with average 
filling n = 0.96875. The dashed line is a power-law decay ob¬ 
tained with an exponent /r = ~Kp and Kp from the Var[lJ] 
and e extrapolations in Table I; the vertical offset is chosen 
to show the agreement with the long-distance behavior of the 
correlation function. 



Figure 12. (color online) Decay of the pair correlation func¬ 
tion D{r) for many system system sizes with average hlling 
n = 0.96875. The dashed line is a power-law decay obtained 
with an exponent v = —IjKp and Kp from the Var[jy] and 
e extrapolations in Table I; the vertical offset is chosen to 
show the agreement with the long-distance behavior of the 
correlation function. 





